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Lewis B. Schiff and Joseph L. Steger+ 
Ames Research Center 


I. SUMMARY 


A noniterative, implicit, space-marching, finite-difference algorithm is 
developed for the steady thin-layer Navier-Stokes equations in conservation- 
law form. The numerical algorithm is applicable to steady supersonic viscous 
flow over bodies of arbitrary shape. In addition, the same code can be used 
to compute supersonic inviscid flow or three-dimensional boundary layers. 
Computed results from two-dimensional and three-dimensional versions of the 
numerical algorithm are in good agreement with those obtained from more costly 
time-marching techniques. 


II. INTRODUCTION 


Considerable effort is being directed toward developing efficient finite- 
difference, numerical algorithms for the solution of the unsteady compressible 
Navier-Stokes equations. Although current algorithms are considerably more 
efficient than those available even a few years ago, the cost of time-marched 
Navier-Stokes solutions is not trivial. Furthermore, the computation of 
viscous flow about practical three-dimensional configurations is currently 
restricted by the size of available computer storage. 

For steady, supersonic, high Reynolds number viscous flows about config- 
urations with moderate axial-geometry variation, a substantial additional 
reduction in both computational effort and required storage can be achieved by 
utilizing the parabolized Navier-Stokes equations. The parabolized Navier- 
Stokes equations are obtained by (1) neglecting the unsteady terms as well as 
the streamwise viscous diffusion terms within the Navier-Stokes equations, and 
(2) by modifying the streamwise convective flux vector to permit stable time- 
like marching of the equations downstream from initial data. The resulting 
equations are commonly referred to as the parabolized Navier-Stokes equations 
because they are parabolic-like with respect to the downstream marching coor- 
dinate. Computational efficiency and reduced storage requirements are 
obtained because the parabolized equations are solved by advancing an initial 
plane of data in space, rather than by advancing an initial cube of data in 
time, as is done for the full Navier-Stokes equations. The parabolized 
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Navier-Stokes approximation has been employed by numerous researchers for both 
external flows (cf. refs. 1-6) and internal flows (cf. refs. 7-12). A variety 
of numerical algorithms has been used to advance the resulting equations and, 
as discussed in section III, many of these are unstable. 

In this report, we present a noniterative. Implicit, finite-difference 
algorithm, analogous to that developed by Beam and Warming (cf. refs. 13, 14) 
for unsteady flow, for the solution of the parabolized Navier-Stokes equations 
The algorithm is conservative, second-order accurate in the marching direction 
and can be second- or fourth-order accurate in the transverse directions. 
Although it was developed Independently, our numerical algorithm is computa- 
tionally similar to one recently reported by Vigneron, Rakich, and Tannehlll 
(ref. 6), but it differs crucially in the treatment of the streamwise pressure 
gradient within the subsonic viscous layer. 

In section III, we detail the governing equations, the parabolized 
Navier-Stokes approximation, and the numerical algorithm for steady two- 
dimensional flow. Section III also contains sample computed results and a 
discussion of the stability of the present method. In section IV, we outline 
the extension of the method to steady three-dimensional flow, and present 
sample results that demonstrate the accuracy and versatility of the resulting 
factored algorithm. 


III. TWO-DIMENSIONAL FLOW 


Discussion of the parabolized Navier-Stokes approximation and illustra- 
tion of the numerical algorithm are facilitated if we first consider the case 
of steady two-dimensional flow. The extension to steady three-dimensional 
flow is given in section IV. 


Governing Equations 

Geneval'lzed coovdi-nate transfoimat-Lon- To accommodate body-conforming 
coordinates, we introduce new independent spatial variables that transform the 
physical x,z plane surrounding the body (fig. 1(a)) into a rectangular 5,? 
computational plane (fig. 1(b)). The transformation, of the form 


5 = C(x) = streamwise coordinate 
C = 4(x,z) = normal coordinate 


( 1 ) 


maps the body surface onto ^ = 0. This transformation both simplifies the 
application of surface boundary conditions, and makes possible the approxima- 
tion of neglecting streamwise viscous terms in high Reynolds number flow. 

The Jacobian of the transformation is 


J 


’x^z 


1 

x^z^ 


( 2 ) 
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(a) Physical plane. 
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Figure 1.- Mapping of physical plane into computational plane. 


and the metric derivatives 5x> ?x> etc., in the computational plane, are 
related to those in the physical plane, x^, z^, etc., by 


In this report we consider 5 to be ?(x) only. Thus, vertical lines in the 
physical plane map into vertical lines of the computational plane. 


Transformation of gasdynamia equations- The steady Navier-Stokes equa- 
tions, written in strong conservation-law form for Cartesian spatial variables 
x,z, can be expressed in nondimensional variables as 


where 


E = E(q) = 


9E 3G 
3x 3z 



= Ji_ + IS \ 

Re \3x 3z/ 


G = G(q) = 


pw 
puw 

pw^ + p 
(e + p)wi 


q = 



(4) 


and the form of the viscous term is discussed in the section on Viscous 
model (p. 4). 


3 


The strong conservation- law form can be preserved under the trans- 

formation of coordinates by retaining the Cartesian velocity components as 
dependent variables (cf. refs. 15, 16), and the transformed equation becomes 


9? 9C Re \95 9?/ (5) 

where 



and the contravariant velocity components U,W are defined in terms of the 
Cartesian velocity components u,w as 


U = ?^u 

W = 


( 6 ) 


The internal energy of the gas, e^, is defined in terms of the conservative 
variables as 


e^ = - 0.5 (u 2 + w2) (7) 

and for a perfect gas with ratio of specific heats y, the equation of state 
is 


P 


P 


(y - l)ej- 


a 


2 


Y 


In equations (4) -(7) the Cartesian velocity components u,w are made non- 
dimensional with respect to aa> (the free-stream speed of sound) , density p 
is normalized by Pcoj and total energy e is referenced to Poo^oo^* 

Viscous model- The first step in obtaining the parabolized Navier-Stokes 
equations from equation (5) is to neglect all streamwise derivatives, 9/95, 
within the viscous terms. This approximation is physically justified for high 
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Reynolds number flow past body-conforming coordinates by using the usual argu- 
ments of boundary-layer theory (see also ref. 17 or 18 for related hypersonic 
viscous-flow analysis) . Neglect of the streamwise viscous terms is necessary 
to prevent exponential growth in marching the equations in that is, to 
mathematically change the nature of equation (5) from elliptic to parabolic 
type with respect to the 5 coordinate. On neglecting the streamwise deriva- 
tives, equation (5) can be written as 

M = J_ li (81 

9? 9? Re 9C ^ ' 


where, in equation (8) 


/ (v/3)(<;xU^ + \ 

S = 7 I + (y/3)(c^u^ + | (9) 

I ?z^) [ (y/2) + w2)^ + KPr"l(Y - l)“^(a2)^] j 

\ + ^^w)(C^u^ + / 

In obtaining equation (9), use has been made of the Stokes hypothesis, 

A = -2y/3, where A and y are the coefficients of viscosity. Also, k is the 
coefficient of thermal conductivity. Re is the free-stream Reynolds number, 
and Pr is the free-stream Prandtl number. For turbulent flow computations 
the eddy-viscosity model described by Baldwin and Lomax (ref. 19) is employed. 

Boundary cond-ltions- Surface boundary conditions for equation (8) are 
simplified because the body surface has been mapped onto ? = 0 (see fig. 1) . 
The steady no-slip condition is simply given by U = W = 0. The pressure on 
the body surface can be determined from the <; momentum equation, evaluated 
at the wall, which becomes 


&■) 


1 _ 9 _ 
Re 9^ 


tx" + (4/3) 


yw^ + 


However, a simplified boundary condition, 9p/9^ =0, is consistent with 
restrictions to be placed on the governing equations to maintain a stable 
streamwise marching procedure. The surface density is obtained from the equa- 
tion of state using the found surface pressure and a specification of either 
the wall temperature or temperature gradient. 

In the present computations no provision has been made for fitting the 
bow shock wave. Instead, the outer edge of the computational region, 

? = Cjnax* chosen to extend into the undisturbed free stream beyond the 
shock layer, and the bow shock is captured. 

A solution consistent with the parabolized Navier-Stokes approximation 
must be supplied as initial data. The initial data must be those of supersonic 
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external flow, and the streamwise component of velocity must be everywhere 
positive. 

Jacob'tan matvi'Oes of the flux vectors- Jacobian matrices of the flux 
vectors are needed in our development of both the parabolized Navier-Stokes 
equations and in the implicit marching algorithm to be described later. Since 
the flux vectors E and G are linear combinations of the Cartesian flux vec- 
tors E and G, 


E 




( 10 ) 


the Jacobian matrices A = [9E/3q] and C = [3G/9q] can be written as 


A = , C = CxA- + ?zC 


( 11 ) 


in terms of the Jacobian matrices of the Cartesian flux vectors A = [3E/9q] 
and C = [3G/3q]. Any one matrix can be obtained from the general form 


K, 


Ko 


0 


A or C = 


Ki(}> 2 - u0 
- w0 

0[2(})2 - Y(e/p)] 


0 - K^(y - 2)u 

Kj^w - K2 (y ~ l)u 

{Kj[Y(e/p) - (f)^] 
- (y - l)u0} 


K2U - (y - l)Kj^w 

0 - K 2 (y - 2)w 
{K2[y(s/p) - (j)2] 
- (y - l)w0} 


Ki(y - 1) 
K2(y - 1) 

Y0 


where 


( 12 ) 


= 0.5(y - 1) (u^ + w^), 0 = K^u + K 2 W 


To obtain the matrix A, set = Cx arid K 2 = 0, and to obtain the Jacobian 
matrix C, set = i^x and K 2 = Cz- The Cartesian Jacobian matrices A and 
C are obtained from equation (12) with = 1, K 2 = 0, and with = 0, 

K 2 = 1, respectively. 


The Cartesian flux vectors E and G are homogeneous functions of degree 
one in q. As a consequence, they possess the identities 

E = Aq G = Cq (13) 

The homogeneous property also extends to the generalized flux vectors, that is, 

E = Aq G = Cq (14) 

With the streamwise variation of the coefficients of viscosity y and 
thermal conductivity k neglected, a Jacobian matrix operator for the viscous 
term can be written as 
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( 15 ) 


with 


0 

0 

0 

0 

™21 

cti3^(l/p) 

a29^(l/p) 

0 

™31 

a23^(l/p) 

a33^(l/p) 

0 

“41 

m42 

“4 3 

“ 49 ^ 


m2i = -ai3^(u/p) - a29^(w/p) 

“31 = -a23^(u/p) - a 38 ^(w/p) 

mm = a43^[-(e/p2) + (u^ + w^)/p] - a^S^Cu^/p) - 2a23^(uw/p) - a33^(w^/p) 
“42 = -a43j-(u/p) - m2i 


“43 = -a43^(w/p) - m3! 

Oil = m[(4/3)Cx^ + » 0^2 = (ij/3)4x4x 

“3 = + ( 4 / 3 )Cz^] » 014 = (YK/Pr)(i;x^ + 4 z^) 


and 

p = (p/J) e = (e/J) 


The viscous term is homogeneous of degree zero in q and thus possesses the 
property that 

Mq = 0 (16) 


The Parabolized Navier-Stokes Approximation 

Cond'Lt'ions for stable marohtng- As alluded to in the introduction, the 
parabolized Navier-Stokes approximation for steady supersonic external flow 
employs two main assumptions: (1) the viscous terms in the marching direction 

5 (which we loosely refer to as streamwise) are negligible, and (2) the stream- 
wise convective flux derivative has positive time-like behavior (discussed 
below) with respect to the remaining spatial derivatives. The first approxi- 
mation is justified for high Reynolds number flow and body-conforming coor- 
dinates and has been discussed above. The second assumption is the most 
difficult restriction to Impose on the parabolized Navier-Stokes equations. 

With the Navier-Stokes equations arranged as in equation (8) , by positive 
time-like behavior in ?, we mean that the Jacobian matrix A has positive 
eigenvalues. Although we oversimplify, the restriction that A has positive 
eigenvalues is required in Inviscid flow for the equations to be hyperbolic 
and it is needed in viscous flow if positive viscosity is to cause damping in 
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the marching direction. Insofar as the viscous flow near a no-slip wall is 
subsonic, at least one eigenvalue of A, the u - a root, will be less than 
zero. Consequently, the solution can grow exponentially with marching unless 
this negative root is suppressed. 


This becomes readily apparent if we consider a linearized frozen (i.e. , 
locally constant) coefficient form of equation (8) . 


Af 


li 

85 


+ Cf 


35 




9^q 

9 + fo 


(17) 


where Af , Cf , and Nf are constant coefficient matrices v/ith elements defined 
by equations (12) and (15). The matrix Nf differs from M insofar as the 
operators 9/35 have been shifted to the right. To the lowest order linear- 
ization = 0, otherwise f^ is a known function and contains linearization 

terms such as Re ^9 ^(Sq - Nf 9^qo) , which result from the expansions 
S = So + Mo(q - qo) j Mq = Nf6^ + e. The metrics are also assumed to be con- 
stant (i.e., uniform grid) as are the coefficients of viscosity and thermal 
conductivity. 


If u =5^ a and u 0, Af^ exists and equation (17) can be rewritten as 


li + = -i_ a-ut 


35 


f 9? Re f 




+ A^"f 


(18) 


The eigenvalues of Af^Cf are (l/5x)(5x "*■ w/u, 

(uw ± a/u^ + w^ - a^)/(u2 - a^)]}. The eigenvalues of Af^Nf, computed with 
the simplification Zx ~ given by 



By introduction of a suitable similarity transform, either matrix product 
Af^Cf or Af^Nf can be diagonalized, but they cannot be simultaneously diag- 
onalized. In particular, if is the matrix whose column vectors are the 

independent eigenvectors of Af^Nf, then introduction of 


X 


-i;:: = 


into equation (18) yields, on premultiplication by X, 

+ So 


35 


+ XAflCfX-1 If = i 


( 20 ) 


( 21 ) 
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where is a diagonal matrix with elements equal to the eigenvalues given 

by equation (19). If 0 < u < a, has one negative real eigenvalue and 

the remaining eigenvalues are positive. This is evident from equation (19) 
where 02 > 0 if u > 0 and the roots cf 3^4 3 ± v^ 3 ^ - 0 (u^ - a^) where 

0 = 4yk/pP^ >0. If u < a, /p + 0(a2 - u^) >3 so one root must be posi- 
tive real and one root must be negative real. Which of the roots 03 or 04 
is negative depends on the magnitude of u. 

Taking 04 as the negative root for u < a, the fourth scalar equation 
of equation (21) is seen to have an effective negative viscosity. As such, 
the fourth vector component ri^ grows exponentially with marching in 
Moreover, even if the diffusion coefficient is negligible, the roots of A^^C^ 
are complex if u^ + w^ < a^ . Thus, the inviscid part is not hyperbolic in 
5 unless u^ + w^ > a^. Consequently, for stable streamwise marching, the 
eigenvalues of must be positive real while the roots of AJ^C^ should 

be real. Although stated as an oversimplification of the matrix algebra, 
these conditions occur precisely when A^ has positive real roots. 

The subsoni-o layer model- Two observations are now made. The first is 
that if pressure can be specified in the flux vector E, that is, the given p 
is not a function of q, then the sound speed contribution to the eigenvalues 
of [3E/9q] = A is removed. In this way, the eigenvalues of A remain posi- 
tive as long as u > 0. The second observation is that, for high Reynolds 
number viscous flow, pressure is approximately constant through the thin sub- 
sonic viscous sublayer near the wall. Indeed, according to boundary-layer 
theory, for high Reynolds number flow the approximation 9p/3n =0 is valid 
over the entire thickness of the viscous layer. Thus, this approximation is 
even more apropos over just the subsonic portion of the viscous layer. 


Although developed under a different formalism, these observations form 
the basis of the parabolized Navler-Stokes approximation that Rubin and Lin 
(cf. refs. 2, 3) term the sublayer approximation. In our development of the 
subsonic layer (i.e., sublayer) approximation, we begin by defining a new 
streamwise flux vector, E^, given by 




( 22 ) 


where Ps = (y - 1) [e - 0.5p(u^ + w^)] for supersonic flow u > a(l + Eg) 
and Pg is defined from 9p/9<; = 0 for subsonic flow u < a(l + £ 3 ). Here 
we assume that c is effectively normal to the surface, and Eg is a small 
positive number picked so that u ^ a and A“^ exists. 


A schematic of how Pg is evaluated is given in figure 2. The essen- 
tial idea is that for points within the subsonic viscous sublayer pg is 
not evaluated from the local flow variables but is taken from the adjacent 




Figure 2.- Schematic of sublayer 
approximation showing pg 
Impressed from above. 


supersonic flow. Throughout the sub- 
sonic sublayer It Is assumed that 
Ps Ps(q)> where q Is the local vec- 
tor of dependent variables. Of course, 
Pg Is related to q In the adjacent 
supersonic region. 


The Jacobian matrix, = [3Eg/3q] 
where pg Is specified, has positive 
real roots If u > 0. By repeating 
the frozen coefficient analysis. It Is 
shown below that equation (8) should be 
stable for marching In 5 if E Is 
replaced by Eg. The Jacobian matrix 
Is given by 


0 


1 0 0 



-uw 


2u 

w 


0 0 
u 0 


(23) 


-u(e + Pg)/p 


(e+ps)/p 0 u 


and has eigenvalues o(Au) = u,u,u,u. Indeed, the vector Eg was originally 
constructed from similarity transforms so that has the eigenvalues of u 

(see ref. 20 for related work). 


The eigenvalues of 


A^^Nf (Au = 5x^u) with Cx = that Is 

(o. 1. I , 


(24) 


are positive real If u > 0. Consequently, according to the frozen coeffi- 
cient analysis, the viscous part of the Initial value problem Is stable for 
marching In We now examine the Invlscld part of the equations.^ If for 

analysis purposes we apply the subsonic sublayer approximation to G, then 
Cy Is defined and all four eigenvalues of A^^Cy are real, are Identical, 
and are given by (l/?x)t4x + Cz(w/u)]. If the sublayer approximation Is not 
Imposed on G, then under the restriction ?x “ we find the eigenvalues of 
A^kzC are 




^ z 

X 


fur w w ± a 1 
Lu ’ u ’ u I 


(25) 


In either case the Invlscld flow has real eigenvalues; hox^ever, the Invlscld 
portion of the parabolized Navler-Stokes equations Is not strictly hyperbolic 
because AJ^^C^ and A-u^C are defective by one eigenvector (as Is A^; see, 
e.g., ref. 21 for a definition of hyperbollclty for first-order systems of 
equations). Apparently the Invlscld part still retains a weak hyperbollclty. 
Curiously, the defective matrices Ay and have eigenvalue properties like 
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those of simple commuting matrices, such as o(Au)a(Cu) = o(AuCu) and 
a(Au) + a(Cu) = o(Au + C^) (note that A^q = uq) . 

Summarizing, the parabolized Navler-Stokes equations with the sublayer 
approximation can be expressed as 

3 E ^ ^ 

35 8? Re 3? 

where Eg Is defined by equation (22) , and the equations are stable for 
marching In 5 when u > 0. 

Relation to other work- Since the parabolized Navler-Stokes equations 
have been used extensively, we feel that It Is Important to show how equa- 
tion (26) Is related to past work. Although usually Investigated without the 
eigenvalue formalism. It has long been recognized (cf. refs. 2, 3, 6, 7) that 
the crucial approximation In all parabolized Navler-Stokes schemes Is In the 
treatment of the p^ (l.e. , 3p/35) term In subsonic flow regions. In our 
development we think In terms of the pressure Itself, but specifying pg as 
a function of 5 is equivalent to specifying p^. 

The earliest successful parabolized Navler-Stokes schemes used the 
approximation that pg = 0 In subsonic regions. This method always proved 
to be stable, as Indeed It should be. According to the frozen coefficient 
theory, the marching should be stable when p^ Is specified In subsonic 
regions and u > 0. Setting p^ = 0 In subsonic regions is equivalent to 
setting Pg to a specified constant, namely, the initial value of Pg. 

Various researchers (refs. 1-6) have attempted to retain p^ in the sub- 
sonic regions. The usual idea has been to lag the differencing of p^ so as 
to treat it as a "source term." This, in fact, simply amounts to an explicit 
differencing of p^. According to the frozen coefficient analysis, however, 
retaining p^ will always lead to instability In the limit of refined grid 
spacing. Thxs occurs because the differential equations themselves admit 
exponential growth in the subsonic region unless the functional dependence of 
p^ Is suppressed. Stability analysis by Lubard and Helliwell (ref. 5) seems 
to suggest that explicit differencing of p^ can lead to weaker Instability 
than Implicit differencing of p^. 

Any scheme that retains p^ in the subsonic sublayer relies on numerical 
dissipation to suppress unstable exponential growth. If the numerical scheme 
is consistent with the differential equations, then according to frozen 
coefficient analysis it will always be unstable as the grid spacing is refined. 

A clever means of using numerical dissipation to control unstable growth 
due to retaining p^ has been employed by Lubard and Helliwell (ref. 5). As 
we interpret their technique, they take advantage of the fact that the 
implicit Euler numerical differencing scheme for initial-value problems will 
be stable in regions where the differential equation itself is unstable (see 
ref. 22 or 23 for the numerical stability bounds on the implicit Euler 
scheme). If the chosen step size is sufficiently large in Ax to suppress 
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the "physical" instability, but not so large as to trigger "nonlinear" insta- 
bility, the Lubard and Helliwell scheme can be used to compute solutions. 
However, the method is always inconsistent in the sense that the grid spacing 
cannot be arbitrarily refined. 

If the Lubard and Helliwell approach seems impractical, it must be 
remarked that the sublayer method also exhibits erratic divergence as the grid 
spacing in Ax is refined. So-called departure solutions are discussed in 
the literature (refs. 2, 3, 6). This behavior and a practical means of con- 
trol will be discussed in the section on Departure Solutions and Global 
Iteration (p. 20). We simply comment here that the sublayer analysis pre- 
viously discussed is a local analysis that does not account for global inter- 
action between Pg and q of the outer flow. 

More recently Vigneron, Raklch, and Tannehill (ref. 6) developed a para- 
bolized Navier-Stokes scheme similar to the one developed in this paper. The 
crucial difference is that in reference 6 the authors attempted to approximate 
the p^ term with a weighting between implicit and explicit differencing that 
depends on the local Mach number. 

Finally, we should note that a variety of parabolized Navier-Stokes 
schemes has been advanced for subsonic internal flow (refs. 7-12). The 
Patankar and Spalding (ref. 7) method appears to be the forerunner of these 
techniques and, as in the sublayer method, Patankar and Spalding determined 
p^ from special contrived relations based on the known mass flux through the 
channel . 


Development of the Numerical Algorithm 


A fully implicit, noniterative, finite-difference algorithm is constructed 
for the parabolized Navier-Stokes equations with the sublayer approximation. 

The difference equations are treated in vector form and their solution requires 
a block tridiagonal inversion at each marching step. Figure 1 indicates the 
extent of the computational domain and the definition of the indices j and 1. 


Difference operators- An implicit, finite-difference scheme for equa- 
tion (26) is constructed by selecting difference operators that would be 
stable for a model problem of diffusion and convection. The following dif- 
ference approximations are selected for the inviscid flux vectors 



(E^^^ - V) - ol(eJ - ^) 


(1 - 


+ 0(A5) 


l+Sct 


where a 


0 or 1/3 for first- or second-order accuracy, and 


9 G . , g ^£- 1-1 ^£-1 

8? ■ ■ 2A? 


0(AC)2 


(27) 


(28) 
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Each term of the viscous flux vector, 8^S, is of the form and is 

differenced as 




^h+1 + ■ V - Vl^ 


2A?2 

Applying these operators to equation (26) gives 


(29) 


(E^'*'^- V) - a(Ej - P b 

(1 - a)A5 


^ - 0 


(30) 


This choice of difference operators is unconditionally stable for the model 
initial-value problem in 5 


3u ^ ^u _ 3^u 
35 “ 3C ~ 9?2 


(31) 


Because A or A^ has positive real roots, we expect the difference equations 
represented by equation (30) to be unconditionally linear stable. 


Loaat tinearizat-ions- To avoid solving a nonlinear system of equations at 
each step in 5, the flux vectors of equation (30) at j+1 are replaced by 
local linearizations about j. The local linearizations are defined as 


,j+l 

^ E^ + Aj(qj + 1 - qj) = A^q^^^ 

(32a) 

J+1 

^ + C^(qj+^ - q^) = 

(32b) 

,:i+i 

- (4+^Aj + cf 'cj)qj + ' E 

(33a) 

J+1 

) 

:: - 4) = 

(33b) 


where all of the approximations are 0(A5)^ and we have used the homogeneous 
property, equation (13). The Jacobian matrices A, C, and M are defined by 
equations (12) and (15). Note that ~ indicates that the matrices C3 , M3 , 
and the vector S3 are evaluated using variables q at j and metric quan- 
tities at j+1 (cf. eq. (33a)). 


The special flux vector Es, u < a, has the functional form Es = Es(q,Ps) 
and locally linearizes as 




J+1 


Ps^) = ‘ + Ps 


J+LJ 


(34) 


where A^ is previously defined _(eq, (23)) while is the vector 

(0,l,0,u)l-. The quantity Pg^^ is also unknown, so we extrapolate 

p^^ = Ps^ + P(Ps^ - Pr^) + 0(^5)^"^^ (35) 
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where g = 0 or 1 for first- or second-order accurate extrapolation. We 
remark that we have only exercised the g = 0 option because the error of the 
first-order approximation, confined to the thin sublayer, has not been signif- 
icant in our test calculations. 


The linearization of 
thus be expressed as 



for either supersonic or subsonic flow can 



(36) 


where 



= a 3 = ri+1 


E ^ = 0 
P 


in supersonic regions, u > a, and 







in subsonic regions, u < a. 

If only E^”P^ is locally linearized, the three-point backward difference 
operator becomes first order and nonconservative. That is, for a = 1/3 


C s 


[Ag-q 






^ + 0(AC)2 - e ^] - a(E ^ 


- eJ-*) 

s 


(1 - a)A? 


(37) 


is 0(A5)^/[(1 - a)(A5)] = O(A^) and is nonconservative. However, if Eg^ 
is also linearized, the lowest-order linearization error will be subtracted 
off, that is. 




+ (5 /J)^’^^E ^ + 0(AC)2 - 

X P s 


i + 1 
S s 


- u /J)^e^ ^ 

X p 


- 0(A^)2] - C((E^^ - E^ b 


(1 - a)A? 


(38) 


The difference approximation, equation (38), is 0(A?)^ and is a conservative 
operator. Even if only first-order accuracy is required in 5, a_= 0 above, 
it is still necessary to linearize each term of the difference - E^J 

to maintain a conservative differencing for shock-capturing purposes. 

Delta form algorithm- Applying the local linearizations and adding a 
fourth-order dissipation term to equation (30) results in 
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) + (1 - a)A5(6^C^ - Re 
s s ^ c, 

= (1 - a)A?Re~l6^sj + a(Egj - E^"^) - [ (5x/J)^"^^Ep^ - (5x/J)^Ej~^] + 

(39) 

The fourth-order dissipation term has the form 


®q^ = e^Ag^ (J (V^A^)^(Jq)^ 


(40) 


and is added to suppress high-frequency oscillations. Here 

= ^l+z ~ ^^l+i + + ^ 1-2 


(41) 


If a = 0, linear-stability analysis of the dissipation term alone indicates 
that Cq must be less than 1/8. 


Finally, the difference equations are put into delta form by subtracting 
[Ag3 + (1 - a)A?(6j'C^ - Re~^6^M3)]q3 from both sides of the equations. The 
finished form of the numerical differencing algorithm is then 



+ (1 - a)A?(6^C^ 
= -(Ag^-Af' 


- Re“l(S^S^} 


-Re-lfi^Mj)](qj + ^- q^) 

)q^ + a(Eg^ - E^"^) - (1 - a) AS{ 6 ^ (E/J) ^ 


(G/J)^] 


(42) 


where we have used _ the fact that M^qJ = 0, while on the right-hand side of 
equation (42), ClqJ, defined by equation (33a), was written in terms of the 
flux vectors. The delta form, in which the left side operates on 
qJ+1 - qj = AqJ , is not as efficient as the nondelta version of the differ- 
ence equation, equation (39). However, the delta form is more convenient in 
three dimensions, and, as discussed below, higher order spatial accuracy is 
easily obtained with the delta form algorithm. 


Note that the flux vector G is not redefined to employ pg in the sub- 
sonic sublayer because experience shows that no inconsistency develops in 
using the conventional definition of G. We remark that q, not E, was used 
throughout as the dependent variable chiefly because it is awkward to express 
G = CA~^E in terms of the special sublayer flux vector. Eg. Real gas 
effects and the viscous terms, especially the turbulent viscosity coefficients, 
are also more conveniently calculated in terms of q. 


For the delta form algorithm, equation (42) , in high Reynolds number 
flow, it is easy to obtain a scheme that is consistent with fourth-order 
accuracy in the Z direction. One simply replaces the second-order right- 
hand side operator 6^ of equation (42) with the conventional five-point 
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LAMINAR 

□ MARCHING CODE 

TIME-DEPENDENT CODE, REF. 13 

TURBULENT 

O MARCHING CODE 
TIME-DEPENDENT CODE, REF. 24 



Figure 5.- Flat plate viscous layer profiles; = 2.0, Rex = 0.832x10®. 

The turbulent marching-code results agree well with the solution obtained by 
the Steger code (ref. 24) using the same turbulence model. 

tlonlifting b'ioonvex a-ivfo-it- The capacity of the marching code to handle 
streamwise variations of geometry was demonstrated by computing the flow over 
a nonlifting, 10% thick, parabolic arc airfoil. As was done for the flat 
plate, a turbulent flow was computed at = 2.0 and Re^ = 1.85x10® 

based on chord c using the time-dependent code (ref. 24). The computational 
grid used is shown in figure 3. Flow-field profiles taken at x/c = 0.10, 

Re^ = 0.185x10® were used as initial data for the marching code, and a 
marching solution was obtained for 0.10 < x/c < 1.0. 

The marching and time-dependent surface-pressure distributions, shown in 
figure 6, demonstrate excellent agreement over the entire airfoil surface. 
Velocity and density profiles through the viscous layer at x/c = 0.90, 

Rejj = 1.6/xlO®, are shown in figure 7. Again, good agreement is observed 
between the marching and time-dependent results. This is not unexpected 
since, as has been discussed, the normal direction spatial-difference 
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operators and the turbulence model of the marching code have the same form as 
those of the time-dependent code. In addition, the time-dependent results 
demonstrate that, at each streamwise station, 9p/3? = 0 through the subsonic 
part of the viscous layer, physically justifying the validity of the viscous 
sublayer approximation made to permit marching. 


Departure Solutions and Global Iteration 

Although the viscous sublayer method has proved to be accurate and versa- 
tile, experience with the sublayer approximation shows that if one continues 
to refine the marching step size. Ax, the method will ultimately diverge. 

This is particularly true unless the initial data are very consistent with the 
sublayer marching equations. The precise cause of the divergence, often 
called a departure solution, is not settled. An intriguing analysis by Lin 
and Rubin (ref. 3) suggests that disturbances can amplify when certain inte- 
gral quantities across the subsonic layer are negative, but it is not clear, 
at least to us, that their analysis sufficiently models the process of impres- 
sing pg from the stable supersonic region. In any event, we find that the 
departure-solution behavior can be controlled by using a global-iteration 
process . 

In the global Iteration technique one initially specifies an entire Pg 
distribution. The sublayer marching method (with Pg specified) is then used 
as part of a relaxation procedure to predict a new flow-field solution and a 
new Pg distribution. The new Pg distribution is then used to obtain an 
improved solution, and so on until the new wall shear stress equals that of 
the previous iteration. Because pg is specified, any small value of Ax 
can be used in the marching scheme. 

A good Initial guess for Pg can be obtained by running the usual sub- 
layer marching procedure with a sufficiently large value of Ax to be stable. 
Alternatively, a constant value of Pg, corresponding to 3pg/3^ = 0, can be 
safely used as an initial guess. Experience with the global iteration tech- 
nique shows that (1) the solution obtained with the viscous sublayer method 
for stable values of Ax are usually accurate, and (2) that even if a poor 
estimate of Pg is initially specified, the global-iteration process is 
rapidly convergent. In most cases, the pressure distribution is converged 
after two iterations, and the skin-friction distribution no longer varies 
after three or four iterations. 

The following computations illustrate the additional stability gained 
from the global-iteration process. The surface-pressure distribution for 
laminar viscous flow over a 10%-thick biconvex airfoil, at = 2.0 and 

Re^o = 1.0x10® based on chord, is shown in figure 8. This solution was 
computed using the time-dependent code of reference 24. The corresponding 
velocity profile through the viscous layer at x/c = 0.8 is shown in figure 9. 
A marching solution, obtained using the sublayer approximation for Pg and 
using Ax = 0.020, is also shown in figures 8 and 9. The results are in 
excellent agreement with the time-dependent solution. The marching solution 
obtained for Ax = 0.010 is identical to that for Ax = 0.020, but when 
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Figure 8.- Surface-pressure distribution on parabolic arc airfoil; Moo ' 

Re^ = 1.0xl0®/c (laminar). 



Figure 9.- Viscous layer velocity profile on parabolic arc airfoil, Moo 

Re„ = 0.8x10^ (laminar), x/c = 0.8. 



Ax = 0.005 was attempted, the solution diverged. However, by using the 
global iteration procedure, with Pg initially constant, rapidly convergent 
solutions (identical to the one obtained for Ax = 0.010) are obtained for the 
smaller step sizes. Ax = 0.005 and Ax = 0.00125. The convergence sequence for 
the wall shear stress distribution over the first three iterations of the 
Ax = 0.00125 case is shown in figure 10. 
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If the time-dependent solution 
was not available, the global itera- 
tion procedure could have been used to 
show that the stable step size 
Ax = 0.020 is sufficiently small to 
maintain accuracy with the sublayer 
method. We remark that as the global 
iteration process is continued in the 
previous example, the wall-shear 
stress distribution remains converged 
to four-place accuracy for the next 
three global iterations. However, 
with continued iteration an oscilla- 
tion will form near the initial pro- 
file unless sufficient spatial damping 
is used. No underrelaxation has been 
used in the Iteration process. 

Our preferred solution technique 
is to use the sublayer approximation 
and not use the global iteration 
scheme. However, if in some part of 
the flow field it becomes necessary to 
refine Ax to check accuracy, the 
above global-iteration process can be 
used economically in that isolated 
region. By employing the iteration 
technique only for isolated segments, 
storage requirements for the Pg dis- 
tributions can be kept negligible. 
Moreover, Pg need not be stored at 
every point in 5 if one is willing 
to use interpolation. 


Boundary-Layer and Invlscid Flow 


Figure 10.- Velocity gradient distri- An added feature of the algorithm 

bution on parabolic arc airfoil developed for the sublayer form of the 

illustrating global iteration parabolized Navier-Stokes equations in 

procedure. conservation-law form is that the same 

computer code can also be used for 

three-dimensional boundary-layer flow (see also ref. 3) or for supersonic 
inviscid flow. One simply has to alter the boundary condition routine and, 
in the case of inviscid flow, not call the viscous subroutines. 
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The surface boundary condition for Inviscid flow is the tangency condi- 
tion, W = 0. Pressure along the body surface can be determined by the normal 
momentum relation 


(5x?x) 'll + + ?z^) = -pU^Cx If + If^ 

All other body-surface quantities needed for the numerical algorithm can be 
obtained by simple extrapolations. Because the equations are in strong 
conservation- law form, shock waves can be captured. However, strong outer bow 
shocks should be fit in hypersonic flow as numerical oscillations near the 
shock wave will likely result in negative pressure and density. All of our 
test calculations have been for supersonic inviscid flow with Moo S 2. As in 
the case of the parabolized Navier-Stokes equations, inviscid supersonic 
marching is only valid about bodies with moderate streamwise variation. The 
calculation should be terminated if a body protuberance (e.g., a canopy) 
generates a significant embedded subsonic region. 

The surface boundary conditions in boundary layer flow are identical to 
those of the parabolized Navier-Stokes equations. For the direct problem, 
that is, p = Pg is specified, all of the necessary boundary-layer edge condi- 
tions can be obtained from the outer inviscid solution with the exception of 
Wg. A relation for Wg is obtained by evaluating the continuity equation at 
the edge, that is 



?=?e 


_3_ 

3 ? 


PeUe 


(45) 


IV. THREE-DIMENSIONAL FLOW 


The development of the implicit marching algorithm for steady three- 
dimensional flow closely parallels the one presented above for two-dimensional 
flow. The same physical assumptions are made, specifically, neglecting the 
streamwise derivatives within the viscous terms, and using the sublayer approx- 
imation. As we shall demonstrate, the Inclusion of the additional spatial 
coordinate leads to a factored sequence of block-tridiagonal equations, whose 
block coefficients are now 5x5 matrices. 

In this section we merely outline the development of the three-dimensional 
algorithm, because the extension from two dimensions is straightforward. 


Transformed Governing Equations 


The three-dimensional steady Navier-Stokes equations, written in nondimen- 
sional form, are 


35 3n 3? Re \35 3n 35/ 


( 46 ) 
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where 


5 = 5(x) = streamwise coordinate 

n = n(x,y,z) = spanwlse (circumferential) coordinate 
C = ^(x,y,z) = normal coordinate 

and the body is assumed to be mapped onto the ? = 0 plane (see fig. 11) . As 
before, we neglect the streamwise derivatives within the viscous terms of 



Figure 11.- Transformation of physical space into computational space. 

equation (46) . This approximation is physically valid for high Reynolds num- 
ber flows, where streamwise-f low gradients within the subsonic viscous layer 
are negligible in comparison with those in the normal direction. The same 
argument permits us also to neglect viscous derivatives along the body in the 
circumferential direction. 

Although it is not necessary to drop the circumferential viscous terms in 
the development of the parabolized Navler-Stokes approximation, doing so sim- 
plifies the computations and is therefore incorporated in the present work. 

The remaining viscous terms, containing only normal derivatives, constitute 
the thin-layer model (cf. refs. 19,^24, 25 for further discussion), and can 
be written as (1/Re) 3S/ 3 i; , where S is given by equation (49) below. Intro- 
ducing the sublayer approximation, the resulting three-dimensional parabolized 
Navier-Stokes equations can be written as 

9^5 , 3F ,3G 1 3S 
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The inviscld flux vectors In equation (47) are 



The thin-layer model viscous term is 

0 

+ (y/3)(^^u^ + i^yV^ + 4^ 

y(42 + cj + C|)v^ + (y/3)(c^u^ + c^v^ + 4^w^) 4y 

S=j| y(42 + 4f + 4i)w^ + (y/3)(?^u^ + 4yV^ + 

{(?x + ?y + s|)f(y/2)(u2 + v2 + w2)j- + KPr“^(y - l)“^(a^)j,] 

+ (y/3)(?jjU + ?yV + ?2w)(?jjUj, 4yVj- + 42^?^ 


( 48 ) 


(49) 


and 


Ps 


= (y “ l)[e “ 0.5p(u^ + v^ + w^)] 


u > a(l + So) 


or 


t=“ 


U < a(l + Eg) 


The metric terms are obtained from chain-rule expansion of x^ , y^, etc. 


and are solved for Py, etc., to give 


= 


=x 

X = 

= J(xrz.) 


n 


,y 

Pz = 


and 




^X = ^(y^^n ■ 

4y = 


(50) 


-= x^(y^z^ - y^z^) 


(51) 


25 



Here the contravariant velocities U, V, and W assume the form 




U = ?^u 

1 






V = + 

riyV + n^w 

► 


(52) 



w = ?^u + 

?yV + J 




The Jacobian matrices 

A, B, and C, 

needed in the 

linearization 

of E, F, 

G, can be written as 







0 

Kl 

K 2 


K 3 

0 


- u 0 

6 - Kj(y- 2 )u 

K 2 U- (y- 1)Kj 

V 

KjU - (y - l)KjW 

Ki(y- 1 ) 

, or C = 

K24|2 - v6 

KiV - K2(y - l)u 

0 - K2(y- 2 )v 


KjV - (y - 1 )K 2 W 

K 2 (y- 1 ) 


Kglj)^ - w 6 

KjW - Kj (y - l)u 

K 2 W- K3<y- 1 ) 

V 

0 - K 3 (y- 2 )w 

K 3 (y- 1 ) 


e[2(j>2 - y(e/p)] 

{Kj [Y(e/p) - 4>^] 

{KjEyCe/p) - ij> 

2] 

{K2[Y(e/p) - <})^] 

Y0 



- (y - l)u 6 } 

- (y-1)v9} 


- (y - l)w 0 1 



(53) 


where = 0.5(y ~ 1) w^) , 0 = Kj^u + K 2 V + K 3 W and, for example, 

to obtain C, Ki = ^2 ~ ?y >' ^3 = Cz- 

The viscous vector S is linearized by Taylor series as in reference 24, 
producing the coefficient matrix operator 


0 

0 

0 

0 

0 


™21 

ai 6 jj(l/p) 

a 2 d^ ( 1 /p) 

a 3 d^(l/p) 

0 


“31 

a 2 d^ ( 1 /p) 

a46^(l/p) 

asd^ ( 1 /p) 

0 

(54) 

“41 

a36jj(l/p) 

asdj- ( 1 /p) 

agdj- ( 1 /p) 

0 


“51 

“52 

“5 3 

“54 

Oio 6 i^(l/p)_ 



with 


“21 

= a^d^(-u/p) 

+ 

a 2 ^l^(-v/p) 

+ 

ot 3 d^(-w/p) 


= U 2 dj,(-u/p) 

+ 

«4'S^(-v/p) 

+ 

a^d^ (-w/p) 

“41 

= agdj.(-u/p) 

+ 

a56^(-v/p) 

+ 

agd^ (-w/p) 

“51 

= a^d^(-u 2 /p) 

+ ^ 2 ^^ (~ 2 uv/p) 

+ a 3 d^(- 2 uw/p) 


+ a^^d^C-v^/p) + ag6^(-w^/p) + agd^(-2vw/p) 

+ agd^C-e/p^) + ag 6 ^[(u^ + v^ + w^)/p] 
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‘52 

= -m2i - aQSj.(u/p) , 

™5 3 

“™31 “ 

‘54 

= -mi^j - ap6^(w/p) 



“o 

= YkPr-l(?2 + ^2 + j.2) ^ 

“1 = 

y[(4/3)c^+ ?|) ] 

“2 

= (y/3)Cx?y . 

“3 = 

(y/3)C2c^z 

“4 

= + ?z^ J > 

“5 = 

(y/3)Cy?2 . 

“6 

= y[4+?2+ (4/3)?|] 




Finally, the sublayer Jacobian matrix i: [3Eg/3q] is given by 

0 




-u‘ 


-uv 


-uw 


1 

2u 

V 


0 0 

0 0 

u 0 

0 u 


0 

0 

0 

0 


-u(e + p^)/p (e + p^)/p 0 


0 


and again all the eigenvalues of are equal to u. 


Numerical Algorithm and Solution Procedure 

The implicit marching algorithm for the solution of equation (47) is 
derived in the same manner as its two-dimensional counterpart, with the 
flux vector linearized in the same manner as G. The resulting algorithm 
written in delta form, is 

[A^^ + (1- a)A?(<S^B^ + - Re" - q^ ) 

S T| tp (a 

= -(Ag^ - +a(Eg^ - E^“^) 

- (1- a)Acj6^[nJ‘^^(E/J)j +n^'^^(F/J)^ +n^‘^^(G/J)j] 

+ 6^[4^^(E/J)^ + + ?J‘^^(G/J)^] - Re-l^^gJj 


(55) 


(56) 
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where 6^ is central differenced like 6^, equation (28), and the smoothing 
term ® is defined by 






Here 


:g must be less than 1/16 for stability. 


An approximately-factored form of equation (56), which retains the same 
order of accuracy in can be obtained if we note that 

+ (1 - a)A5(6^B^)j (Ag^)“^[Ag^ + (1 - a)AC^6^C^ - Re“^6^M^)j Aq^ 


= LHS(56) + 0(A5)3 (57) 


(Note: Ag^ can degrade the factorization error if u is sufficiently small.) 

On replacing the left side of equation (56), LHS(56), with the left side of 
equation (57), one obtains the factored algorithm. The algorithm is solved by 
the sequence of implicit inversions 


pgJ + (1 - a)AC(6^BJ)]Aq* = RHS(56) 


+ (1 - a)A^^6^C^ - Re"^6^M^^jAq^ = 


As^Aq, 


(58a) 

(58b) 


Equation (58) differs from its two-dimensional analogy primarily in the inclu- 
sion of the implicit circumferential inversion factor. 


A typical computational grid is shown in the physical crossflow plane 
[x = Xq(Cq)] in figure 12. The grid extends radially between the body 

surface and an outer boundary 


k = 10 


k = KMAX 


k= 19 


e = LMAX 



Figure 12.- Cross-section of typical 
computational grid, x = Xq(5q)* 


located in the undisturbed free 
stream, and is chosen to completely 
circumscribe the body, 1 ^ k ^ KMAX, 
to permit a treatment of 
nonbilaterally-symmetric flows. 

Such flows include the case of com- 
bined angles of attack and yaw, 
and the important case of a non- 
symmetric leeward side wake exhib- 
k = 28 ited by axisymmetrlc bodies at 

large incidence. S 3 nnmetric flows 
can be treated with half the com- 
putational effort by employing a 
grid that runs from the windward to 
leeward plane of symmetry and by 
applying the usual symmetry condi- 
tions at the edges. 


28 


To advance the solution of equation (58) , we first form the right-hand 
side terms of equation (58a) and perform the circumferential Implicit inver- 
sion- The use of a central difference approximation for the n derivatives, 
together with the periodic continuation condition, leads to a periodic block- 
trldiagonal system of equations. This system is inverted, using the solver 
described in reference 26, to obtain the intermediate variables. Once these 
quantities are known, the right-hand side terms of equation (58b), AgdAq^, 
are evaluated, and the equation is inverted in the normal direction, using the 
same procedure previously described for equation (42) , to obtain AqJ . , and 
thus qj'*’^ . ’ 


Three-Dimensional Results 

The accuracy of the factored marching algorithm applied to three- 
dimensional flow was evaluated by computing the flow field about a hemisphere- 
cylinder body at 0° and at 5° angle of attack. The test-case conditions were 
again chosen to duplicate steady flow-field results obtained from time- 
dependent Navier-Stokes computations and, for the body at incidence, to match 
those of the wind-tunnel experiment described in reference 27. 

Ax'Lsyrnmetvic f'Low- Although the flow field surrounding the hemisphere 
cylinder at zero incidence is axisymmetric , the Cartesian velocity components 
used in the computation vary sinusoidally in the circumferential direction 
around the body. Thus, this case provides a nontrivial test of the factoriza- 
tion procedure. The azimuthal-invariant time-dependent code, described in 
reference 28, was used to compute the turbulent flow around the body, at 
= 2.0 and Re^o = 8.80x10^ based on nose radius R]j;j, using the grid 
shown in longitudinal section in figure 13. The flow field exhibits an 
embedded subsonic region in the shock 
layer at the nose, which expands around 
the nose and becomes supersonic in the 
vicinity of the sphere-cylinder junc- 
tion. Flow profiles taken at 
x/Rj^ = 3.45, downstream of the subsonic 
region, were used as initial data for 
the marching code, and a marching solu- 
tion was obtained for 3.45 < x/Rj^ < 21.0 

The marching and time-dependent 
surface-pressure distributions are shown 
in figure 14 and are in good agreement 
over the entire body. The small axial 
oscillation in the marching results is 
attributed to a small inconsistency 
between the initial data and the march- 
ing technique. The amplitude of the 
oscillation is never more than 1% of the 
pressure and is seen to damp toward the 
rear of the body. Velocity and density 
profiles within the viscous layer. 



x/R|yi 


Figure 13.- Axisymmetric hemisphere 
cylinder computational mesh. 


29 


1.1 



0 2 4 6 8 10 12 14 16 18 20 22 

x/R|yj 

Figure 14.- Axisynmetrlc hemisphere cylinder surface-pressure distribution; 

= 2.0, Re^ = 8 . SOxlO^^/Rjg (turbulent). 

taken from the marching and time-dependent solutions at x/R^^j = 20.6, are 
shown in figure 15. The marching results show good agreement with those of 
the time-dependent code. 

Hemisphere cylinder at incidence- The flow field surrounding a hemi- 
sphere cylinder at incidence in a low-Mach-number supersonic stream has 
recently been investigated experimentally by Hsieh (ref. 27), and computa- 
tionally by Pulliam and Steger (ref. 25), who used a three-dimensional, time- 
dependent, thin-layer Navier-Stokes code. Their computational grid was 
selected to resolve the details of the flow in the region of the nose, and in 
this region the computed results are in good agreement with the experimental 
measurements. However, the limitation of computer storage required that the 
grid be progressively stretched axially along the cylinder. Consequently, the 
streamwise details of the downstream flow were only marginally resolved. The 
use of the marching code, with initial data taken from the time-dependent 
solution in a region of good resolution, can circumvent the storage limita- 
tions. Using a grid similar to that shown in figure 13, a steady turbulent 
flow solution was obtained using the time-dependent code (ref. 25) at 
= 1.40, Re^ = 2.0x10^ based on Rj>j, and a = 5°. Data taken at 
x/Rj;f = 3.07 were prescribed as initial data and a marching solution was 
obtained from 3.07 < x/Rj^ < 40.0. A comparison of the surface-pressure dis- 
tributions along the windward and leeward planes of symmetry is shown in fig- 
ure 16, together with the experiments of Hsieh (ref. 27). Although the 
marching solution was obtained for x/Rj^ < 40.0, and could be continued down- 
stream, only the data for the region where the marching results, the time- 
dependent results, and the experimental measurements overlap, 3 . 07 < x/Rj^ < 16 . 0, 
are presented in figure 16. The marching results are in good agreement with 
the time-dependent results in the region common to both computations, 
x/R^ < 14.0, and both are in good agreement with the measured surface pres- 
sures. However, the marching-code results give better agreement with the 
measured values for 9.0 < x/Rj^ < 14.0, where the time-dependent solution 
lacks resolution. 
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u/Uoo p/Pcx= 


Figure 15.- Axisymmetric hemisphere cylinder viscous layer profiles; Moo = 2.0, 

Rex = 1.81x10^ (turbulent), x/Rj.j = 20.6. 

Streamwise velocity profiles through the viscous layer on the windward 
and leeward rays, taken from the computational results at x/Rjj = 6.98, are 
shown in figure 17. At this axial location, the stretched grid of the time- 
dependent solution still maintains adequate streamwise resolution. The veloc- 
ity gradient is much more sensitive than is surface pressure. Thus, the good 
agreement between the time-dependent and the marching solutions attests to the 
accuracy of the factored marching algorithm. Also, the time-dependent results 
exhibit constant pressure across the subsonic viscous layer, thus justifying 
the assumptions made in the viscous sublayer approximation. 


Conical Flow Fields 

Motivation- The initial data for the marching method must, in general, be 
supplied from an auxiliary, time-dependent computation. However, when consid- 
ering flow over conical or pointed bodies, the marching technique can be used 
to generate its own Initial data. For example, inviscld conical solutions can 
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Figure 16.- Windward and leeward symmetry plane surface-pressure distributions 
on hemisphere cylinder at incidence; Moo = 1.40, Re^ = 2.0xl0^/Rj;[ (turbu- 
lent) , a = 5 ° . 



Figure 17.- Viscous layer velocity profiles on hemisphere cylinder at incidence; 
Moo = 1.40, Re^j. = 1.40x10^ (turbulent), x/Rj^ = 6.98, a = 5°. 
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be obtained by using the marching method as a distance-asymptotic technique. 
The computational grid is chosen to be conical, with grid points at successive 
axial stations located along rays emanating from the cone apex (see fig. 18). 
The flow variables are initially set 
to free-stream values and the equa- 
tions are marched downstream from 
X = Xq to X = Xq + Ax. After each 
step, the solution is scaled to place 
it back at the initial station, 

X = Xq. When no change in the flow 
variables occurs with further march- 
ing, the flow variables are constant 
along rays, and a conical flow field 
has been obtained. 

If conditions within the viscous 
layer are also assumed to be conical 
(see ref. 29 for discussion and 
ref. 30 for experimental confirmation 
for high Reynolds number flows) , the 
marching step-back procedure can be 
used to generate conical viscous 
flows. Here, the assumption of coni- 
cal flow permits setting 9pg/3£; = 0 
within the subsonic viscous layer. 

In this case, the marching step-back 
method is numerically equivalent to 
that of McRae (ref. 29), in which the 
Navier-Stokes equations are written 
in conical coordinates, derivatives 
along rays are dropped, and the 
resulting equations are advanced in 
time to obtain a steady solution. 

Con-ieaZ results- A series of computations was performed to obtain laminar 
flows over a 9.09° half-angle cone at = 2,0, Re^ = 1,85x10^ based on 
axial distance from the nose, and for angles of attack ranging from 0° to 15°. 
The computational grid completely encircled the body and the resulting flows 
were found to be symmetric about the windward plane. 

The circumferential surface-pressure distribution found for a = 10° is 
shown in figure 19, and is in good agreement with the corresponding results 
obtained by McRae (private communication, AFFDL, Ames Research Center, Moffett 
Field, Calif., 1978). Pressure contours in the crossflow plane, x = Xq, of 
the marching solution are shown in figure 20, and demonstrate the symmetry of 
the flow. 

At 10° angle of attack, a small reversed crossflow-separation region 
occurs near the leeward plane of symmetry. This can be seen in figure 21, 
which presents the projections of the flow velocity vectors onto the crossflow 
plane for points near the body surface. The location of the circumferential 



Figure 18.- Conical flow grid. 
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Figure 19.- Circumferential surface-pressure distribution on 9.09° half-angle 
cone; Moo = 2.0, Re^ = 1.85x10^ (laminar), a = 10°. 



0 .018 Y 


Figure 20.- Crossflow plane pressure Figure 21.- Crossflow plane velocity 
contours; 0^ = 9.09°, Moo = 2.0, vectors; 0£, = 9.09°, Moo = 2.0, 

Rex = 1.85x10^ (laminar), a = 10°. Rex = 1.85x10^ (laminar), a = 10°. 


separation point 0g, obtained by interpolation, is also indicated in 
figure 21. 

The circumferential surface-pressure distribution and crossflow velocity 
vectors for a = 15° are shown in figures 22 and 23, respectively. At this 
angle of attack, relative incidence a/d^ = 1.65, the crossflow separation 
region is more pronounced radially. Also the crossflow separation point is 
located closer to the windward symmetry plane. 



Figure 22.- Circumferential surface-pressure distribution; 0^ = 9.09°, 
Mco = 2.0, Re^ = 1.85x10^ (laminar), a = 15°. 
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Figure 23.- Crossflow plane velocity vectors; 6^ = 9.09°, Moo = 2.0, 

Re^^ = 1.85x10^ (laminar), a = 15°. 


V. CONCLUDING REMARKS 


A noniterative, implicit, finite-difference marching algorithm has been 
developed for steady supersonic viscous flow. The parabolized Navier-Stokes 
equations, in strong conservation-law form, have been transformed into general 
coordinates so that arbitrary body shapes can be mapped onto constant planes 
in the uniform computational space. The approximately factored finite- 
difference algorithm is noniterative, second-order accurate in the marching 
direction, and second- or fourth-order accurate in the crossflow plane. Use 
of the subsonic layer approximation with a global iteration technique for 
"surface" pressure allows the grid spacing to be refined in a uniform manner. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif. 94035, May 8, 1980 
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